remove(list=ls())



#Load necessary packages: reader, lme4, stragazer, ggplot2, ggeffects

library(readr)
library(lme4)
library(stargazer)
library(ggeffects)
library(MuMIn)


#Load the data "burdensharing_dat.csv"



burdensharing_dat<-read.csv("burdensharing_dat.csv")


#############################################
#Models reported in Table 2 in the main text#
#############################################

#Model 1

m1lmer<-lmer(trps_prop_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy + cnl + other_mission + mountains +  (1|mission_id) + (1|contribccode), dat=subset(burdensharing_dat, osv_2years==1 | conflict_2years==1))
summary(m1lmer)

r.squaredGLMM(m1lmer)

#Model 2

m2lmer<-lmer(trps_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy  + cnl + other_mission + mountains +  (1|mission_id) + (1|contribccode), dat=subset(burdensharing_dat, osv_2years==1 | conflict_2years==1))
summary(m2lmer)

r.squaredGLMM(m2lmer)



stargazer(m1lmer, m2lmer, title="Results from Multi-level Linear Regression", 
          align=T, digits=4, intercept.bottom=T, dep.var.labels=c("Proportion of Troops", "Count of Troops"), star.cutoffs=c(0.1, 0.05, 0.01, 0.001), no.space=T, type="text")




#Marginal effects 

ggpredict(m1lmer, "ln_refugees[0, 6]")

ggpredict(m2lmer, "ln_refugees[0, 6]")

ggpredict(m1lmer, "milex_pers_hth[0.01, 0.79]")

ggpredict(m2lmer, "milex_pers_hth[0.01, 0.79]")

ggpredict(m2lmer, "democracy[0, 1]")


#########################################################
#Summary Statistics Reported in Table 1 in the main text#
#########################################################

stargazer(m1lmer@frame, digits=4, summary.stat=c("mean", "median", "sd", "min", "max"), style="io", type="text", report="vcp")
stargazer(m2lmer@frame, digits=4, summary.stat=c("mean", "median", "sd", "min", "max"), style="io", type="text", report="vcp")



#############################################
#Models reported in Table 3 in the Appendix #
#############################################

#Model 3 

m3lm_fixedfx<-lm(trps_prop_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy + cnl +  as.factor(mission_id), dat=subset(burdensharing_dat, osv_2years==1 | conflict_2years==1))

summary(m3lm_fixedfx)

#Model 4

m4lm_fixedfx<-lm(trps_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy + cnl +  as.factor(mission_id), dat=subset(burdensharing_dat, osv_2years==1 | conflict_2years==1))

summary(m4lm_fixedfx)

stargazer(m3lm_fixedfx, m4lm_fixedfx, 
          align=T, digits=4, intercept.bottom=T, dep.var.labels="Count of Troops", star.cutoffs=c(0.1, 0.05, 0.01, 0.001), no.space=T, type="text")

#############################################
#Models reported in Table 4 in the Appendix #
#############################################


#Model 5

m5lmer_pre2000<-lmer(trps_prop_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy + cnl + other_mission + mountains + (1|mission_id) + (1|contribccode), dat=subset(burdensharing_dat, year<2000 & (osv_2years==1 | conflict_2years==1)))
summary(m5lmer_pre2000)

r.squaredGLMM(m5lmer_pre2000)

#Model 6

m6lmer_post1998<-lmer(trps_prop_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy + cnl + other_mission + mountains +  (1|mission_id) + (1|contribccode), dat=subset(burdensharing_dat, year>1998 & (osv_2years==1 | conflict_2years==1)))
summary(m6lmer_post1998)

r.squaredGLMM(m6lmer_post1998)

#Model 7

m7lmer_pre2000<-lmer(trps_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy  + cnl + other_mission + mountains +  (1|mission_id) + (1|contribccode), dat=subset(burdensharing_dat, year<2000 & (osv_2years==1 | conflict_2years==1)))
summary(m7lmer_pre2000)

r.squaredGLMM(m7lmer_pre2000)

#Model 8

m8lmer_post1998<-lmer(trps_confosv ~ contig_dum +  ln_refugees + milex_pers_hth + democracy + cnl + other_mission + mountains +  (1|mission_id) + (1|contribccode), dat=subset(burdensharing_dat,  year>1998 & (osv_2years==1 | conflict_2years==1)))
summary(m8lmer_post1998)

r.squaredGLMM(m8lmer_post1998)


stargazer(m5lmer_pre2000, m6lmer_post1998, m7lmer_pre2000, m8lmer_post1998, 
          align=T, digits=4, intercept.bottom=T, dep.var.labels="Count of Troops", star.cutoffs=c(0.1, 0.05, 0.01, 0.001), no.space=T, type="text")